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A power spectral density estimation method and 
apparatus . 

TECHNICAL FIELD 



The present invention relates to a bias compensated spectral 
estimation method and apparatus based on a parametric auto- 
regressive model . 

BACKGROUND OF THE INVENTION 

The present invention may be applied, for example, to noise 
suppression [1, 2] in telephony systems, conventional as well as 
cellular, where adaptive algorithms are used in order to model 
and enhance noisy speech based on a single microphone measure- 
ment . 

Speech enhancement by spectral subtraction relies on, explicitly 
or implicitly, accurate power spectral density estimates 
calculated from the noisy speech. The classical method for 
obtaining such estimates is periodograra based on the Fast Fourier 
Transform (FFT) . However, lately another approach has been 
suggested, namely parametric power spectral density estimation, 
which gives a less distorted speech output, a better reduction of 
the noise level and remaining noise without annoying artifacts 
("musical noise") . For details on parametric power spectral 
density estimation in general, see [3, 4]. 

In general, due to model errors, there appears some bias in the 
spectral valleys of the parametric power spectral density 
estimate. In the output from a spectral subtraction based noise 
canceler this bias gives rise to an undesirable "level pumping" 
in the background noise. 

SUMMARY OF THE INVENTION 



An object of the present invention is a method and apparatus that 
eliminates or reduces this "level pumping" of the background 
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noise with relatively low complexity and without numerical 
stability problems. 

This object is achieved by a method and apparatus in accordance 
with the enclosed claims. 



5 The key idea of this invention is to use a data dependent (or 

adaptive) dynamic range expansion for the parametric spectrum 
model in order to improve the audible speech quality in a 
spectral subtraction based noise canceler. 



BRIEF DESCRIPTION OF THE DRAWINGS 



The invention, together with further objects and advantages 
thereof , may best be understood by making reference to the 
following description taken together with the accoTi5)anying 
drawings, in which: 



is a block diagram illustrating an embodiment of an 
apparatus in accordance with the present invention ,- 

is a block diagram of another embodiment of an 
apparatus in accordance with the present invention,- 

is a diagram illustrating the true power spectral 
density, a parametric estimate of the true power 
spectral density and a bias compensated estimate of 
the true power spectral density; 



is another diagram illustrating the true power 
spectral density, a parametric estimate of the true 
power spectral density and a bias compensated 
estimate of the true power spectral density; 



FIGURE 5 



is a flow chart illustrating the method performed 
by the embodiment of Fig. 1; and 
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FIGURE 6 is a flow chart illustrating the method performed 

by the embodiment of Fig . 2 . 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

Throughout the drawings the same reference designations will be 
used for corresponding or similar elements. 

Furthermore, in order to simplify the description of the present 
invention, the mathematical background of the present invention 
has been transferred to the enclosed appendix. In the following 
description numerals within parentheses will refer to correspon-, 
ding equations in this appendix. 

Figure 1 shows a block diagram of an embodiment of the apparatus 
in accordance with the present invention. A frame of speech 
{x(k) } is forwarded to a LPC analyzer (LPC analysis is described 
in, for exan^jle, [5] ) . LPC analyzer 10 determines a set of filter 
coefficients (LPC parameters) that are forwarded to a PSD 
estimator 12 and an inverse filter 14 . PSD estimator 12 determi- 
nes a parametric power spectral density estimate of the input 
frame {x(k)} from the LPC parameters (see (1) in the appendix). 
In Fig. 1 the variance of the input signal is not used as an 
input to PSD estimator 12. Instead a unit signal "l" is forwarded 
to PSD estimator 12. The reason for this is simply that this 
variance would only scale the PSD estimate, and since this 
scaling factor has to be canceled in the final result (se (9) in 
the appendix) , it is simpler to eliminate it from the PSD 
calculation. The estimate from PSD estimator 12 will contain the 
"level pumping" bias mentioned above. 

In order to compensate for the "level pumping" bias the input 
frame {x(k) } is also forwarded to inverse filter 14 for forming 
a residual signal (see (7) in the appendix) , which is forwarded 
to another LPC analyzer 16. LPC analyzer 16 analyses the residual 
signal and forwards corresponding LPC parameters (variance and 
filter coefficients) to a residual PSD estimator 18, which forms 
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a parametric power spectral density estimate of the residual 
signal (see (8) in the appendix) . 

Finally the two parametric power spectral density estimates of 
the input signal and residual signal, respectively, are multi- 
plied by each other in a multiplier 2 0 for obtaining a bias 
compensated parametric power spectral density estimate of input 
signal frame {x(k)} (this corresponds to equation (9) in the 
appendix) . 

Example 

The following scenario is considered: The frame length N=1024 and 
the AR (AR=AutoRegressive) model order p=10 . The underlying true 
system is modeled by the ARMA (ARMA=AutoRegressive- Moving 
Average) process 

x{k) = 1-^-^^0.9^-^ _ ^J^^ 

1-3 .64Z-2-4 .44^-^ ■•-2 . 62z-*-0 .77 2"= 

where eCk) is white noise. 



Figure 3 shows the true power spectral density of the above 
process (solid line) , the biased power spectral density estimate 
from PSD estimator 12 (dash-dotted line) and the bias compensated 
power spectral density estimate in accordance with the present 
invention (dashed line) . From Fig. 3 it is clear that the bias 
compensated power spectral density estimate in general is closer 
to the underlying true power spectral density. Especially in the 
deep valleys (for example for w/ (2Tr) <=0 . 17) the bias compensated 
estimate is much closer (by 5 dB) to the true power spectral 
density. 

In a preferred embodiment of the present invention a design 
parameter y may be used to multiply the bias compensated 
estimate. In Fig. 3 parameter y was assumed to be equal to 1. 
Generally -y is a positive number near i . In the preferred 
embodiment 7 has the value indicated in the algorithm section of 
the appendix. Thus, in this case y differs from frame to frame. 
Pig. 4 is a diagram similar to the diagram in Fig. 3, in which 
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the bias compensated estimate has been scaled by this value of -y . 

, The above described embodiment of Fig. l may be characterized as 

a frequency domain compensation, since the actual compensation is 
performed in the frequency domain by multiplying two power 
5 spectral density estimates with each other. However, such an 

operation corresponds to convolution in the time domain. Thus, 
there is an equivalent time domain iirplementation of the 
invention. Such an embodiment is shown in Fig. 2. 

In Pig. 2 the input signal frame is forwarded to LPC analyzer 10 

10 as in Fig. 1. However, no power spectral density estimation is 

performed with the obtained LPC parameters. Instead the filter 
parameters from LPC analysis of the input signal and residual 
signal are forwarded to a convolution circuit 22, which forwards 
the convoluted parameters to a PSD estimator 12 ' , which forms the 

15 bias compensated estimate, which may be multiplied by 7. The 

convolution step may be viewed as a polynomial multiplication, in 
which a polynomial defined by the filter parameters of the input 
signal is multiplied by the polynomial defined by the filter 
parameters of the residual signal. The coefficients of the 

20 resulting polynomial represent the bias compensated LPC-parame- 

ters. The polynomial multiplication will result in a polynomial 
of higher order, that is, in more coefficients. However, this is 
no problem, since it is customary to "zero pad" the input to a 
PSD estimator to obtain a sufficient number of samples of the PSD 

25 estimate. The result of the higher degree of the polynomial 

obtained by the convolution will only be fewer zeroes. 

Flow charts corresponding to the embodiments of Figs. 1 and 2 are 
given in Figs. 5 and S, respectively. Furthermore, the correspon- 
ding frequency and time domain algorithms are given in the 
► 30 appendix. 

A rough estimation of the numerical complexity may be obtained as 
follows. The residual filtering (7) requires -Np operations (sum 
+ add) . The LPC analysis of e(k) requires «Np operations to form 
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the covariance elements and <=p^ operations to solve the corre- 
sponding set of equations (3) . Of the algorithms (frequency and 
time domain) the time domain algorithm is the most efficient, 
since it requires ^p^ operation for performing the convolution. 
To summarize, the bias compensation can be performed in =2p(N+p) 
operations/frame. For example, with n=256 and p=iO and 50% frame 
overlap, the bias compensation algorithm requires approximately 
0,5x10* instructions/s . 

In this specification the invention has been described with 
reference to speech signals. However, the same idea is also 
applicable in other applications that rely on parametric spectral 
estimation of measured signals. Such applications can be found, 
for example, in the areas of radar and sonar, economics, optical 
interferometry, biomedicine, vibration analysis, image pro- 
cessing, radio astronomy, oceanography, etc. 

It will be understood by those skilled in the art that various 
modifications and changes may be made to the present invention 
without departure from the spirit and scope thereof, which is 
defined by the appended claims . 
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APPENDIX 

Clonsidcr the real-valued zero mean signal {x(A;)}, k = 1 N where N denotes the 

frame length {N = 160, for example). The autoregressive spectral estimator (ARSPE) is 
given by, see |3, 4| 



where u is the angular frequency uj £ (0, 27r). In (1), A(z) is given by 



(1) 



(2) 



where 0^ = (o i • - - Op)^ are the estimated AR coefficients (found by LPC analysis, see 
|5|) and is the residual error variance. The estimated parameter vector and are 
calculated from {.t(A:)} as follows: 

4 = -R-'f 



= ^ x{t + k)x(i) 



(5) 



The set of linear equations (3) can be solved using the Levinson-Durbin algorithm, see 
13]. The spectral estimate (1) is known to be smooth and its statistical properties have 
been analyzed in [6] for broad-band and noisy narrow-band signals, respectively. 

In general, due to model errors there appears some bias in the spectral valleys. Roughly, 
this bias can be described as 

I « 0 for w such that ct>i(cj) w max^ 
> 0 for w such that «: max<„ ^^(f^) 

where is the estimate (1) and 4»i(i<^) is the true (and unknown) power spectral 

density oi x{k). 



(6) 
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In order r,o reduce the bias appearing in the spectral valleys, the residual is calculated 



according to 

e{k) = A{z-')x{k) 



Performing another LPC analysis on {£•(/:)}, the residual power spectral density can be 
ciilculated from. cf. (1) 

where, similarly to (2), 0. = (6i • ■ -fig)"^ denotes the estimated AR coefficients and a; the 
error variance. In general,the model order q p, but here it seems reasonable to let p = g. 
Preferably p ~ \//V, for example N may be chosen around 10. 

In the proposed frequency domain algorithm below, the estimate (1) is compensated 
according to 

4>x{u;) = ^ • «l>e(c^) • <f.,(a;) (9) 

where 7 {w 1 ) is a design variable. The frequency domain algorithm is summarized in the 
algorithms section below and in the block diagrams in Fig. 1 and 5. 

A corresponding time domain algorithm is also summarized in the algorithms section and 
in Fig. 2 and 6. In this case the compensation is performed in a convolution step, in 
which the LPC filter coefficients Ox are compensated. This embodiment is more efficient, 
since one PSD estimation is replaced by a less complex convolution. In this embodiment 
the scaUng factor 7 may simply be set to a constant near or equal to 1. However, it 
is also possible to calculate 7 for each frame, as in the frequency domain algorithm by 
calculating the root of the characteristic polynomial defined by 9^ that lies closest to the 
unit circle. If the angle of this root is denoted uJ , then 
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ALGORITHMS 

INPUTS 

X input data X = 
p LPC model order 

OUTPUTS 

Or signal LPC parametei-s 0:r. = (oi • • • a^V 

(Tj signal LPC residual variance 

4>.^ signal LPC spectrum = (<J>i(l ) ■ ■ ■ i'^iN/2)f 

$3, compensated LPC spectrum 4>x = • • - 4>i(Ar/2))^ 

£ residual c = (e(l)---f(Ar))^ 

§c residual LPC parameters 9^ = (6i - • • 6^)^ 

(?j residual LPC error variance 

7 design variable (=l/{max)t ^^(fc)) in preferred embodiment) 
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FREQUENCY DOMAIN ALGORITHM 

FOR EACH FRAME DO THE FOLLOWING STEPS: 
(power spectra! density estimation ) 
\0x, &l\ ■■= LPCanalyze(x, p] signal LPC analysis 
^:r. SPECC^j:^, \, N) Signal spectral estimation, a'i set to 1 

(bias compensation) 
e ■- FILTER(^:„ x) residual filtering 

|4, := LPCanalyze( e, p) residual LPC analysis 
4>e := SPEC(^E, a;, N) residual spectral estimation 

FOR k=l TO NI2 DO spectral compensation 

^^{k) := 7 ■ ^r(k) ■ i>e(^-) l/(max*,*^(jk)) < 7 < 1 
END FOR 

TIME DOMAIN ALGORITHM 

FOR EACH FRAME DO THE FOLLOWING STEPS: 

\Qx, al] := LPCanalyze(x, p) signal LPC analysis 

e := FILTER(^^, x) residual filtering 

l^j, ajl := LPCanalyze(£, p) residual LPC analysis 

e ■.=CONV{9^,9^) LPC compensation 

^ := SPEC(e", a;, N) spectral estimation 
FOR k=l TO N/2 DO 

*i(fc) := 7 • |)(fc) scaling 
END FOR 
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12 
CLAIMS 

1. A power spectral density estimation method, coitiprising the 
steps of : 

performing a LPC analysis on an input signal vector for 
determining a first set of LPC filter parameters; 

determining a first power spectral density estimate of said 
input signal vector based on said first set of LPC filter parame- 
ters ; 

filtering said input signal vector through an inverse LPC 
filter determined by said first set of LPC filter parameters for 
obtaining a residual signal vector,- 

performing a LPC analysis on said residual signal vector for 
determining a second set of LPC filter parameters; 

determining a second power spectral density estimate of said 
residual signal vector based on said second set of LPC filter 
parameters ; and 

forming a bias compensated power spectral estimate of said 
input signal vector that is proportional to the product of said 
first and second power spectral estimates. 

2. The method of claim 1, wherein said product is multiplied by 
a positive scaling factor that is less than or equal to 1. 

3. The method of claim 2, wherein said scaling factor is the 
inverted value of the maximum value of said second power spectral 
density estimate. 

4. The method of claim 1, 2 or 3, wherein said input signal 
vector comprises speech samples. 

5. A power spectral density estimation method, comprising the 
steps of : 

performing a LPC analysis on an input signal vector for 
determining a first set of LPC filter parameters; 

filtering said input signal vector through an inverse LPC 
filter determined by said first set of LPC filter parameters for 
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obtaining a residual signal vector; 

performing a LPC analysis on said residual signal vector for 
, determining a second set of LPC filter parameters,- 

convolving said first set of LPC filter parameters with said 
^ 5 second set of LPC filter parameters for forming a compensated set 

of LPC filter parameters; 

determining a bias compensated power spectral density estimate 
of said input signal vector based on said compensated set of LPC 
filter parameters. 

10 6 . The method of claim 5 , wherein said bias compensated power 

spectral density estimate is multiplied by a positive scaling 
factor that is less than or equal to l . 

7. The method of claim 6, wherein said scaling factor is the 
inverted value of the maximum value of a power spectral density 

15 estimate of said residual signal vector. 

8. The method of claim 5, 6 or 7, wherein said input signal 
vector comprises speech saitples . 

9. A power spectral density estimation apparatus, comprising: 
means (10) for performing a LPC analysis on an input signal 

vector for determining a first set of LPC parameters; 

means (12) for determining a first power spectral density 
estimate of said input signal vector based on said first set of 
LPC parameters; 

means (14) for filtering said input signal vector through an 
inverse LPC filter determined by said first set of LPC parameters 
for obtaining a residual signal vector; 

means (16) for performing a LPC analysis on said residual 
signal vector for determining a second set of LPC parameters; 

means (18) for determining a second power spectral density 
estimate of said residual signal vector based on said second set 
of LPC parameters; and 

means (20) for forming a bias compensated power spectral 
estimate of said input signal vector that is proportional to the 



25 
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product of said first and second power spectral estimates. 

10. A power spectral density estimation apparatus, comprising: 
means (10) for performing a LPC analysis on an input signal 

vector for determining a first set of LPC filter parameters; 

means (14) for filtering said input signal vector through an 

inverse LPC filter determined by said first set of LPC filter 

parameters for obtaining a residual signal vector; 

means (16) for performing a LPC analysis on said residual 

signal vector for determining a second set of LPC filter 

parameters ; 

means (22) for convolving said first set of LPC filter 
parameters with said second set of LPC filter parameters for 
forming a compensated set of LPC filter parameters; 

means (12' ) for determining a bias cott^jensated power spectral 
density estimate of said input signal vector based on said 
compensated set of LPC filter parameters. 
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